

# AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- 
# AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- AEP.frame -- 

setClass("AEP.frame", representation(
	data = "data.frame",
	pobs_data = "matrix",
	nb_obs = "numeric",
	nb_var = "numeric",
	EmpiricalCDF = "list",
	ECDF = "function",
	N = "numeric",	
	binary_matrix = "array",
	sum_ones = "array",
	Hmultiplier = "array",
	alpha = "numeric",
	nb_iteration = "numeric",
	m = "array"
	)
)

setMethod("initialize", "AEP.frame", function(.Object, data, alpha, iter=1) {
	
  if (missing(data)) stop("syntax error")
	if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'") 	  
  
	if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  
	if (is.numeric(iter) == FALSE) stop("the parameter 'iter' must be of type 'numeric'") 	 
  if (iter < 1) stop("the parameter 'iter' must be >= 1")			
  
	.Object@data <- data
	.Object@nb_obs <- dim(data)[1]
	.Object@nb_var <- dim(data)[2]  
	.Object@alpha <- alpha	
	.Object@nb_iteration <- iter
	.Object@N <- 2^.Object@nb_var
	.Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
	.Object@sum_ones <- array(0, c(.Object@N))
	.Object@Hmultiplier <- array(0, c(.Object@N))
			
	#Compute the ECDF of each variable
	EmpiricalCDF <- NULL	
	for (i in 1:.Object@nb_var) .Object@EmpiricalCDF <- c(.Object@EmpiricalCDF, ecdf(.Object@data[,i]))
	.Object@ECDF <- function(i, x) .Object@EmpiricalCDF[[i]](x)

  # Compute the pobs of data
  .Object@pobs_data <- pobs(.Object@data)
  #.Object@pobs_data <- array(0, c(.Object@nb_obs, .Object@nb_var))
  #for (i in 1:.Object@nb_obs) {
  #  for (j in 1:.Object@nb_var) {
  #    .Object@pobs_data[i,j] <- .Object@ECDF(j, .Object@data[i,j]) 
  #  }      
  #}
  #.Object@pobs_data <- as.matrix(.Object@pobs_data)
  
  #Initialize the matrix of binary and the sum of ones in each line
	for (i in 2:.Object@N) {
		.Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
		j <- 1
		.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
		while (.Object@binary_matrix[i, j] > 1) {
			.Object@binary_matrix[i, j] <- 0
			j <- j + 1
			.Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
		}
		.Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
	}

	#Initialize the multiplier of H()
	for (i in 1:.Object@N) {
		if ((.Object@nb_var  -  .Object@sum_ones[i]) %% 2 == 0) .Object@Hmultiplier[i] <- 1
		else  .Object@Hmultiplier[i] <- -1
	}

	# Initialize the vector m
	.Object@m <- array(0, c(.Object@N))
	.Object@m[1] <- 0
	for (j in 2:.Object@N) {
		if (.Object@sum_ones[j] < 1 / .Object@alpha) {
			if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
			else .Object@m[j] <- -1
		}
		else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
		else {
			if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
			else .Object@m[j] <- -1
		}
	}
	.Object
})

as.aep.frame <- function(data, alpha, iter=1) {
  if (missing(data)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'") 	    
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  if (is.numeric(iter) == FALSE) stop("the parameter 'iter' must be of type 'numeric'") 	 
  if (iter < 1) stop("the parameter 'iter' must be >= 1")	    
  new("AEP.frame", data, alpha, iter)
}

setMethod(f = "show", signature = "AEP.frame", definition = function(object) {
  cat("\n@nb_obs = number of observations = ", object@nb_obs, "\n\n")
  cat("@nb_var = number of variables = ", object@nb_var, "\n\n")
  cat("@data = [not displayed]\n\n")
	#print(object@data)
  cat("@pobs_data (pobs(data)) = [not displayed]\n\n")
  #print(object@pobs_data)
  cat("@ECDF(i, x) = empirical cumulative distribution function of var 'i' at 'x', i.e. P[var_i <= x]\n\n")
  cat("@N = pow(2, nb_var) = ", object@N, "\n\n")
  cat("@binary_matrix = \n")
  print(object@binary_matrix)
  cat("\n@sum_ones = sum of each row of the binary matrix = ", object@sum_ones, "\n\n")
  cat("@Hmultiplier = (-1)^(nb_var - sum_ones) = ", object@Hmultiplier, "\n\n")
  cat("@alpha = ", object@alpha, "\n\n")
  cat("@nb_iteration = ", object@nb_iteration, "\n\n")
  cat("@m =", object@m, "\n")  
})


# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 
# AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- AEP.result -- 

setClass("AEP.result", representation(
  prob = "numeric",
  vsum = "numeric",
  nb_simplex = "numeric"
  )
)
setMethod("initialize", "AEP.result", function(.Object, prob, vsum, nb_simplex) {
  .Object@prob <- prob
  .Object@vsum <- vsum
  .Object@nb_simplex <- nb_simplex  
  .Object
})
setMethod(f = "show", signature = "AEP.result", definition = function(object) {
  cat("P[sum <= ", object@vsum ,"] = ", object@prob, "\n")
  cat("Number of simplex resolved = ", object@nb_simplex, "\n")
})


# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 
# AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- AEP.HCSimplex -- 

setClass("AEP.HCSimplex", representation(
	b = "array",
	h = "numeric",
	s = "numeric",
	iteration = "numeric"
	)
)
setMethod("initialize", "AEP.HCSimplex", function(.Object, or, lim, multiplier, it) {
  .Object@b <- or
	.Object@h <- lim
	.Object@s <- multiplier
	.Object@iteration <- it
	.Object
})


# Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- 
# Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- Functions -- 

VHmeasure <-function(frame, hc) {
	if (missing(frame) || missing(hc)) stop("syntax error")
	if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
	if (class(hc) != "AEP.HCSimplex") stop("the parameter 'hc' must be of type 'AEP.HCSimplex'")
	if (length(hc@b) != frame@nb_var) stop("hc@b is too long")		
	
	samplePoints <- array(0, c(frame@N, frame@nb_var))
	Uval <- array(0, c(frame@N, frame@nb_var))
	for (i in 1:frame@N) {
		samplePoints[i,] <- hc@b + (hc@h * frame@binary_matrix[i,])
		for (j in 1:frame@nb_var) Uval[i,j] <- frame@ECDF(j, samplePoints[i,j])	
	}
	H <- C.n(u=Uval, U=frame@pobs_data)
	Hm <- frame@Hmultiplier * H	
	abs(sum(Hm))
}

resolveAEP <- function(frame, vsum) {
  if (missing(frame) || missing(vsum)) stop("syntax error")
  if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")  
	prob <- 0
  nb_simplex <- 0
	spx <- new("AEP.HCSimplex", array(0, c(frame@nb_var)), vsum, 1, 1)
	ListOfSimplexes <- c(spx)	
	while (length(ListOfSimplexes) > 0) {
		spx <- ListOfSimplexes[[1]]
		if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
		else ListOfSimplexes <- NULL
		hc <- new("AEP.HCSimplex", spx@b, frame@alpha * spx@h, spx@s, spx@iteration)
		if (spx@iteration < frame@nb_iteration) {
			new.it <- spx@iteration + 1
			for (j in 2:frame@N) {
				new.b <- spx@b + (frame@alpha * spx@h * frame@binary_matrix[j,])
				new.h <- (1 - (frame@alpha * frame@sum_ones[j])) * spx@h
				new.s <- spx@s * frame@m[j]				
				ListOfSimplexes <- c(ListOfSimplexes, new("AEP.HCSimplex", new.b, new.h, new.s, new.it))			
			}
		}
		prob <- prob + hc@s * VHmeasure(frame, hc)
		nb_simplex <- nb_simplex + 1
	}
  new("AEP.result", prob, vsum, nb_simplex)
} 

AEP.d <- function(data, alpha, iter=1, vsum) {
	if (missing(data) || missing(vsum)) stop("syntax error")
	if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'")    
	if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")        
	if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
	else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
	if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
	if (is.numeric(iter) == FALSE) stop("the parameter 'iter' must be of type 'numeric'") 	 
	if (iter < 1) stop("the parameter 'iter' must be >= 1")	
	frame <- as.aep.frame(data, alpha, iter) 
	resolveAEP(frame, vsum)
}

AEP.f <- function(frame, vsum) {    
	if (missing(frame) || missing(vsum)) stop("syntax error")
	if (class(frame) != "AEP.frame") stop("the parameter 'frame' must be of type 'AEP.frame'")
	if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")
	resolveAEP(frame, vsum)
}

# -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 
# -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 
# -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 -- AEP2 

setClass("AEP.frame2", representation(
  data = "data.frame",
  pobs_data = "matrix",
  nb_obs = "numeric",
  nb_var = "numeric",
  EmpiricalCDF = "list",
  ECDF = "function",
  N = "numeric",	
  binary_matrix = "array",
  sum_ones = "array",
  Hmultiplier = "array",
  alpha = "numeric",
  nb_simplex = "numeric",
  m = "array"
)
)

setMethod("initialize", "AEP.frame2", function(.Object, data, alpha, nb_simplex=1) {
  
  if (missing(data)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'") 	  
  
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  
  if (is.numeric(nb_simplex) == FALSE) stop("the parameter 'nb_simplex' must be of type 'numeric'") 	 
  if (nb_simplex < 1) stop("the parameter 'nb_simplex' must be >= 1")			
  
  .Object@data <- data
  .Object@nb_obs <- dim(data)[1]
  .Object@nb_var <- dim(data)[2]  
  .Object@alpha <- alpha	
  .Object@nb_simplex <- nb_simplex
  .Object@N <- 2^.Object@nb_var
  .Object@binary_matrix <- array(0, c(.Object@N, .Object@nb_var))			
  .Object@sum_ones <- array(0, c(.Object@N))
  .Object@Hmultiplier <- array(0, c(.Object@N))
  
  #Compute the ECDF of each variable
  EmpiricalCDF <- NULL	
  for (i in 1:.Object@nb_var) .Object@EmpiricalCDF <- c(.Object@EmpiricalCDF, ecdf(.Object@data[,i]))
  .Object@ECDF <- function(i, x) .Object@EmpiricalCDF[[i]](x)
  
  # Compute the pobs of data
  .Object@pobs_data <- pobs(.Object@data)
  #.Object@pobs_data <- array(0, c(.Object@nb_obs, .Object@nb_var))
  #for (i in 1:.Object@nb_obs) {
  #  for (j in 1:.Object@nb_var) {
  #    .Object@pobs_data[i,j] <- .Object@ECDF(j, .Object@data[i,j]) 
  #  }      
  #}
  #.Object@pobs_data <- as.matrix(.Object@pobs_data)
  
  #Initialize the matrix of binary and the sum of ones in each line
  for (i in 2:.Object@N) {
    .Object@binary_matrix[i,] <- .Object@binary_matrix[i-1,]
    j <- 1
    .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    while (.Object@binary_matrix[i, j] > 1) {
      .Object@binary_matrix[i, j] <- 0
      j <- j + 1
      .Object@binary_matrix[i, j] <- .Object@binary_matrix[i, j] + 1
    }
    .Object@sum_ones[i] <- sum(.Object@binary_matrix[i,])
  }
  
  #Initialize the multiplier of H()
  for (i in 1:.Object@N) {
    if ((.Object@nb_var  -  .Object@sum_ones[i]) %% 2 == 0) .Object@Hmultiplier[i] <- 1
    else  .Object@Hmultiplier[i] <- -1
  }
  
  # Initialize the vector m
  .Object@m <- array(0, c(.Object@N))
  for (j in 2:.Object@N) {
    if (.Object@sum_ones[j] < 1 / .Object@alpha) {
      if((1 + .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
    else if (.Object@sum_ones[j] == 1 / .Object@alpha) .Object@m[j] <- 0
    else {
      if((.Object@nb_var + 1 - .Object@sum_ones[j]) %% 2 == 0) .Object@m[j] <- 1
      else .Object@m[j] <- -1
    }
  }
  .Object
})

as.aep.frame2 <- function(data, alpha, nb_simplex=1) {
  if (missing(data)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'") 	    
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  if (is.numeric(nb_simplex) == FALSE) stop("the parameter 'nb_simplex' must be of type 'numeric'") 	 
  if (nb_simplex < 1) stop("the parameter 'nb_simplex' must be >= 1")	    
  new("AEP.frame2", data, alpha, nb_simplex)
}

setMethod(f = "show", signature = "AEP.frame2", definition = function(object) {
  cat("\n@nb_obs = number of observations = ", object@nb_obs, "\n\n")
  cat("@nb_var = number of variables = ", object@nb_var, "\n\n")
  cat("@data = [not displayed]\n\n")
  #print(object@data)
  cat("@pobs_data (pobs(data)) = [not displayed]\n\n")
  #print(object@pobs_data)
  cat("@ECDF(i, x) = empirical cumulative distribution function of var 'i' at 'x', i.e. P[var_i <= x]\n\n")
  cat("@N = pow(2, nb_var) = ", object@N, "\n\n")
  cat("@binary_matrix = \n")
  print(object@binary_matrix)
  cat("\n@sum_ones = sum of each row of the binary matrix = ", object@sum_ones, "\n\n")
  cat("@Hmultiplier = (-1)^(nb_var - sum_ones) = ", object@Hmultiplier, "\n\n")
  cat("@alpha = ", object@alpha, "\n\n")
  cat("@nb_simplex = ", object@nb_simplex, "\n\n")
  cat("@m =", object@m, "\n")  
})

VHmeasure2 <-function(frame, hc) {
  if (missing(frame) || missing(hc)) stop("syntax error")
  if (class(frame) != "AEP.frame2") stop("the parameter 'frame' must be of type 'AEP.frame2'")
  if (class(hc) != "AEP.HCSimplex") stop("the parameter 'hc' must be of type 'AEP.HCSimplex'")
  if (length(hc@b) != frame@nb_var) stop("hc@b is too long")		
  
  samplePoints <- array(0, c(frame@N, frame@nb_var))
  Uval <- array(0, c(frame@N, frame@nb_var))
  for (i in 1:frame@N) {
    samplePoints[i,] <- hc@b + (hc@h * frame@binary_matrix[i,])
    for (j in 1:frame@nb_var) Uval[i,j] <- frame@ECDF(j, samplePoints[i,j])	
  }
  H <- C.n(u=Uval, U=frame@pobs_data)
  Hm <- frame@Hmultiplier * H	
  abs(sum(Hm))
}

resolveAEP2 <- function(frame, vsum) {
	if (missing(frame) || missing(vsum)) stop("syntax error")
	if (class(frame) != "AEP.frame2") stop("the parameter 'frame' must be of type 'AEP.frame2'")
	if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")  
	prob <- 0
	nb_simplex <- 0
	spx <- new("AEP.HCSimplex", array(0, c(frame@nb_var)), vsum, 1, 1)
	ListOfSimplexes <- c(spx)
	#while (length(ListOfSimplexes) > 0 & nb_simplex < frame@nb_simplex - length(ListOfSimplexes)) {
	while (length(ListOfSimplexes) > 0 & nb_simplex < frame@nb_simplex) {
		spx <- ListOfSimplexes[[1]]    
		if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
		else ListOfSimplexes <- NULL
		hc_full <- new("AEP.HCSimplex", spx@b, spx@h, spx@s, spx@iteration)
		prob_full <- VHmeasure2(frame, hc_full)		
		#nb_simplex <- nb_simplex + 1
		if (prob_full > 0) {
			hc <- new("AEP.HCSimplex", spx@b, frame@alpha * spx@h, spx@s, spx@iteration)
			prob <- prob + hc@s * VHmeasure2(frame, hc)			
			nb_simplex <- nb_simplex + 1
			new.it <- spx@iteration + 1
			for (j in 2:frame@N) {
				new.b <- spx@b + (frame@alpha * spx@h * frame@binary_matrix[j,])
				new.h <- (1 - (frame@alpha * frame@sum_ones[j])) * spx@h
				new.s <- spx@s * frame@m[j]				
				ListOfSimplexes <- c(ListOfSimplexes, new("AEP.HCSimplex", new.b, new.h, new.s, new.it))			
			}			
		}    
	}      
	#while (length(ListOfSimplexes) > 0) {
	#	spx <- ListOfSimplexes[[1]]
	#	if (length(ListOfSimplexes) > 1) ListOfSimplexes <- ListOfSimplexes[2:length(ListOfSimplexes)]
	#	else ListOfSimplexes <- NULL
	#	hc_full <- new("AEP.HCSimplex", spx@b, spx@h, spx@s, spx@iteration)
	#	if (hc_full@s > 0) prob <- prob + hc_full@s * VHmeasure2(frame, hc_full)      
	#	nb_simplex <- nb_simplex + 1
	#}
	new("AEP.result", prob, vsum, nb_simplex)
} 

AEP.d2 <- function(data, alpha, nb_simplex=1, vsum) {
  if (missing(data) || missing(vsum)) stop("syntax error")
  if (is.data.frame(data) == FALSE) stop("the parameter 'data' must be of type 'data.frame'")    
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")        
  if (missing(alpha)) alpha <- 2 / (dim(data)[2] + 1)
  else if (is.numeric(alpha) == FALSE) stop("the parameter 'alpha' must be of type 'numeric'")      
  if (alpha < 1 / dim(data)[2] || alpha >= 1) stop("the parameter 'alpha' must be in [1/d, 1)")   		
  if (is.numeric(nb_simplex) == FALSE) stop("the parameter 'nb_simplex' must be of type 'numeric'") 	 
  if (nb_simplex < 1) stop("the parameter 'nb_simplex' must be >= 1")	
  frame <- as.aep.frame2(data, alpha, nb_simplex) 
  resolveAEP2(frame, vsum)
}

AEP.f2 <- function(frame, vsum) {    
  if (missing(frame) || missing(vsum)) stop("syntax error")
  if (class(frame) != "AEP.frame2") stop("the parameter 'frame' must be of type 'AEP.frame2'")
  if (is.numeric(vsum) == FALSE) stop("the parameter 'vsum' must be of type 'numeric'")
  resolveAEP2(frame, vsum)
}

# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 
# TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- TEST -- 

test.debug <- FALSE

if (test.debug) {
  
	#test case
  nb.var <- 2
	nb.obs <- 30000	
	mydata = matrix(runif(nb.obs * nb.var), c(nb.obs, nb.var))
  for (i in 1:15000) {    
    mydata[i,1] <- mydata[i,1] * 5
    mydata[i,2] <- mydata[i,2] * 5 + 20 
  }
  for (i in 15001:30000) {    
    mydata[i,1] <- mydata[i,1] * 5 + 20
    mydata[i,2] <- mydata[i,2] * 5  
  }
	
	mydata <- data.frame(mydata)
  
  plot(mydata, xlim=c(0,25), ylim=c(0,25))
  somme <- mydata[,1] + mydata[,2]
  print(sum(somme <= 25) / length(somme))
  
  nb_iter <- 8 
	#print(AEP.d(data=mydata, vsum=25, iter=nb_iter))    
  test.frame <- as.aep.frame(data=mydata, iter=nb_iter)  
  print(AEP.f(frame=test.frame, vsum=25))
  
  # Computation of the number of simplexes solved by AEP.x
  nb_splx <- 0
  for (i in 1:nb_iter) nb_splx <- nb_splx + (2^(dim(mydata)[2]) - 1)^(i-1)
  print(paste("nb simplexes = ", nb_splx))
  
  # Utilization of AEP.x2 with the same target number of simplexes
  #print(AEP.d2(data=mydata, vsum=25, nb_simplex=nb_splx))
  test.frame <- as.aep.frame2(data=mydata, nb_simplex=nb_splx)  
  print(AEP.f2(frame=test.frame, vsum=25))
  
	#end of test case
  
} else {	
  
  par(ask=FALSE)
	par(mfrow=c(1,1))
	mydata <- read.table("analysis_traces.txt", header=TRUE, sep=" ")
  nb.obs <- dim(mydata)[1]
  nb.var <- dim(mydata)[2]
  
  # mydata cannot contain any 0 !!!!!
  # otherwise the estimation is screwed because of the P[X_k  in (b_k, b_k + h]] in equation 1.3. The problem is the "(" with b_k not included.
  # try with 2 variables X_1 and X_2 and the following values:
  # (10, 0) and (0, 50)
  # For Q( (0,0) , 50 ) we get H(0,0) - H(0, 50) - H(50, 0) + H(50, 50) = 0 - 0.5 - 0.5 + 1 = 0 !
  
  mydata[mydata == 0] <- 0.0001
  
	observed.exectime <- array(0, c(nb.obs))
	for (i in 1:nb.obs) observed.exectime[i] <- sum(mydata[i,])
	plot( x = 1:nb.obs, y = observed.exectime, 
		main = "Observed execution time",
		xlab = "#run",
		ylab = "Observed Exec Time",
		type = "l", 
		col = "red")

	observed.ecdf <- ecdf(observed.exectime)
	plot(observed.ecdf,
		main = "Observed ECDF",
		xlab = "#run",
		ylab = "P[ExecTime < X]"
	)
  
  	
	my.step <- (max(observed.exectime) - min(observed.exectime)) / 50
	estimated.ecdf <- array(0, c(50))

	my.aep.frame <- as.aep.frame2(data=mydata, nb_simplex=32)
  
	index <- 1
	steps <- seq(min(observed.exectime), max(observed.exectime), by=my.step)
	for (i in steps) {
    result <- AEP.f2(frame=my.aep.frame, vsum=i)
    print(result)
		estimated.ecdf[index] <- result@prob
		index <- index + 1
	}
	
  for (i in seq(length(estimated.ecdf), 1, -1)) {
    if (estimated.ecdf[i] > 1) estimated.ecdf[i] <- 1
    if (i < length(estimated.ecdf) & estimated.ecdf[i] > estimated.ecdf[i+1]) estimated.ecdf[i] <- estimated.ecdf[i+1]
  }
  
	lines( x = steps, y = estimated.ecdf, 		
		type = "l", 
		col = "red")
  
  
  
  print("now processing AEP")
  
  my.step <- (max(observed.exectime) - min(observed.exectime)) / 50
  estimated.ecdf <- array(0, c(50))
  
  my.aep.frame <- as.aep.frame(data=mydata, iter=3)

  index <- 1
  steps <- seq(min(observed.exectime), max(observed.exectime), by=my.step)
  for (i in steps) {  
    result <- AEP.f(frame=my.aep.frame, vsum=i)  
    print(result)  
    estimated.ecdf[index] <- result@prob
    index <- index + 1  
  }
    
  for (i in seq(length(estimated.ecdf), 1, -1)) {
    if (estimated.ecdf[i] > 1) estimated.ecdf[i] <- 1
    if (i < length(estimated.ecdf) & estimated.ecdf[i] > estimated.ecdf[i+1]) estimated.ecdf[i] <- estimated.ecdf[i+1]
  }
  
  lines( x = steps, y = estimated.ecdf,          
         type = "l", 
         col = "green")
  
}
